;***************************************************************
;discription： ndvi/npp与气温、降水的偏关系数、复相关系数以及t检验，f检验函数分析程序入点
;根据-京津冀地区NDVI变化及气候因子驱动分析_孟丹-很重要@相关系数公式完成；
;本程序完成 ndvi/npp与气温、降水的偏关系数、复相关系数以及t检验，f检验函数分析运算工作
;
;author: wangqiang cdut 2015-10-29
;email：  cqwensten@163.com  419436974@qq.com
;***************************************************************
pro pro_rela_calculate
  ;严格编译规则
  COMPILE_OPT IDL2
  ;二次开发模式初始化，恢复envi的sav文件
  ENVI, /RESTORE_BASE_SAVE_FILES
  ;初始化envi，保存日志文件
  ENVI_BATCH_INIT,log_file='relativity.log'

  ;是否退出程序
  isexit = 0
  starttime = systime(1)
  PRINT , '-----' ,string(starttime), '-----'
  ;********************************开始****************************************


  npp_file= envi_pickfile(title='选择NDVI/NPP文件');这个函数可进行文件选择
  ;open the image file
  envi_open_file, npp_file, r_fid=npp_fid ; open the ndvi file
  if (npp_fid eq -1) then begin
    isexit = 1
  endif

  tmp_file= envi_pickfile(title='选择气温文件');这个函数可进行文件选择r
  ;open the image file
  envi_open_file, tmp_file, r_fid=tmp_fid ; open the ndvi file
  if (tmp_fid eq -1) then begin
    isexit = 1
  endif

  pre_file= envi_pickfile(title='选择降水文件');这个函数可进行文件选择
  ;open the image file
  envi_open_file, pre_file, r_fid=pre_fid ; open the ndvi file
  if (pre_fid eq -1) then begin
    isexit = 1
  endif

  if (isexit eq 1) then begin
    endtime = systime(1)
    duration = double(endtime)-string(starttime)
    endstr = '-----持续时间:' +string(endtime)+ '-----'
    PRINT ,endstr

    Dlg_Result= Dialog_Message('输入的数据不完整，相关分析退出(the data is not saitifed)！',title='提示',/Information)
    ;销毁资源

    ;关闭envi二次开发模式
    ENVI_BATCH_EXIT
    return
  endif


  ;获取行、列号，波段数，数据类型
  envi_file_query, npp_fid, ns=ncols, nl=nrows, nb=nbands, Data_Type=data_type, xstart=xstart, ystart=ystart
  print,'the npp file info(npp_fid,列数,行数,波段数) :',npp_fid, ncols,nrows, nbands

  npp_data = FltArr(ncols,nrows,nbands)
  tmp_data = FltArr(ncols,nrows,nbands)
  pre_data = FltArr(ncols,nrows,nbands)

  ;获取文件投影
  map_info = envi_get_map_info(fid=npp_fid)
  help,map_info
  print,map_info

  ;为了更好看到所做的工作，打印npp_file输入路径
  print,npp_file

  dims = [-1, 0, ncols - 1, 0, nrows - 1 ] ;roi_fid,列的起始位置，列的终止位置，行的起始位置，行的终止位置
  FOR i=0,nbands-1 DO BEGIN
    ;get the image data in the openning file which band pos=i
    npp= envi_get_data(fid=npp_fid, dims=dims, pos=i) ;波段数从0开始,提取出训练数据,并且第一层数据为目标数据
    npp_data[*,*,i] = npp

    tmp= envi_get_data(fid=tmp_fid, dims=dims, pos=i) ;波段数从0开始,提取出训练数据,并且第一层数据为目标数据
    tmp_data[*,*,i] = tmp

    pre= envi_get_data(fid=pre_fid, dims=dims, pos=i) ;波段数从0开始,提取出训练数据,并且第一层数据为目标数据
    pre_data[*,*,i] = pre
  ENDFOR

  ;1、计算person相关系数
  rel_npptmp = func_relativity_b2b(ncols,nrows,npp_data,tmp_data) ;npp与气温的person相关系数
  rel_npppre = func_relativity_b2b(ncols,nrows,npp_data,pre_data)  ;npp与降水的person相关系数

  ;2、计算偏相关系数（净相关系数）
  prel_npptmp = func_prel_b2b(ncols,nrows,npp_data,tmp_data,pre_data) ;NPP与气温的偏相关系数
  prel_npppre = func_prel_b2b(ncols,nrows,npp_data,pre_data,tmp_data) ;NPP与降水的偏相关系数

  ;3、计算复相关系数
  frel_npptmppre = func_frel_b2b(ncols,nrows,npp_data,tmp_data,pre_data)  ;NPP与气温、降水的复相关系数

  prel_npptmp[where(finite(prel_npptmp,/NAN)) ]  = -999.0 ;移除NaN为-9999.0
  prel_npppre[where(finite(prel_npppre,/NAN)) ]  = -999.0 ;移除NaN为-9999.0
  frel_npptmppre[where(finite(frel_npptmppre,/NAN)) ]  = -999.0 ;移除NaN为-9999.0

  ;保存模型结果
  output_file = DIALOG_PICKFILE(TITLE = 'NDVI/npp与气温的person相关系数...')
  IF strlen(strtrim(output_file,2)) GT 0 THEN BEGIN
    ;ENVI_WRITE_ENVI_FILE, prel_npptmp, OUT_NAME = output_file,Map_Info = map_info
    ;二进制方式输出
    OPENW,lun, output_file,/get_lun
    ;写入npp结果数据到文件中
    WRITEU,lun,rel_npptmp
    FREE_LUN,lun
    ;写出头文件信息
    ENVI_SETUP_HEAD,fname=output_file, ns=ncols, nl=nrows, nb=1, interleave=0, data_type=4, offset=0, map_info=map_info, /write, /open
  ENDIF

  output_file = DIALOG_PICKFILE(TITLE = 'NDVI/npp与降雨的person相关系数...')
  IF strlen(strtrim(output_file,2)) GT 0 THEN BEGIN
    ;二进制方式输出
    OPENW,lun, output_file,/get_lun
    ;写入npp结果数据到文件中
    WRITEU,lun,rel_npppre
    FREE_LUN,lun
    ;写出头文件信息
    ENVI_SETUP_HEAD,fname=output_file, ns=ncols, nl=nrows, nb=1, interleave=0, data_type=4, offset=0, map_info=map_info, /write, /open
  ENDIF

  output_file = DIALOG_PICKFILE(TITLE = 'NDVI/npp与气温的偏相关系数...')
  IF strlen(strtrim(output_file,2)) GT 0 THEN BEGIN
    ;二进制方式输出
    OPENW,lun, output_file,/get_lun
    ;写入npp结果数据到文件中
    WRITEU,lun,prel_npptmp
    FREE_LUN,lun
    ;写出头文件信息
    ENVI_SETUP_HEAD,fname=output_file, ns=ncols, nl=nrows, nb=1, interleave=0, data_type=4, offset=0, map_info=map_info, /write, /open
  ENDIF

  ;保存模型结果
  output_file = DIALOG_PICKFILE(TITLE = 'NDVI/npp与降水的偏相关系数...')
  IF strlen(strtrim(output_file,2)) GT 0 THEN BEGIN
    ;二进制方式输出
    OPENW,lun, output_file,/get_lun
    ;写入npp结果数据到文件中
    WRITEU,lun,prel_npppre
    FREE_LUN,lun
    ;写出头文件信息
    ENVI_SETUP_HEAD,fname=output_file, ns=ncols, nl=nrows, nb=1, interleave=0, data_type=4, offset=0, map_info=map_info, /write, /open
  ENDIF

  ;保存模型结果
  output_file = DIALOG_PICKFILE(TITLE = 'NDVI/npp与气温、降水的复相关系数...')
  IF strlen(strtrim(output_file,2)) GT 0 THEN BEGIN
    ;二进制方式输出
    OPENW,lun, output_file,/get_lun
    ;写入npp结果数据到文件中
    WRITEU,lun,frel_npptmppre
    FREE_LUN,lun
    ;写出头文件信息
    ENVI_SETUP_HEAD,fname=output_file, ns=ncols, nl=nrows, nb=1, interleave=0, data_type=4, offset=0, map_info=map_info, /write, /open
  ENDIF


  alpha_t = sqrt(nbands -1 - 1) ;sqrt(n-m-1)
  tchk_npptmp = func_check_t(ncols, nrows, alpha_t, prel_npptmp)
  tchk_npppre = func_check_t(ncols, nrows, alpha_t, prel_npppre)

  alpha_f = float(nbands -2 -1)/2.0 ;(n-k-1)/k
  fchk_npptmppre = func_check_f(ncols, nrows, alpha_f, frel_npptmppre)


  ;保存模型结果
  output_file = DIALOG_PICKFILE(TITLE = 'NDVI/npp与气温的偏相关系数T检验结果...')
  IF strlen(strtrim(output_file,2)) GT 0 THEN BEGIN
    ;二进制方式输出
    OPENW,lun, output_file,/get_lun
    ;写入npp结果数据到文件中
    WRITEU,lun,tchk_npptmp
    FREE_LUN,lun
    ;写出头文件信息
    ENVI_SETUP_HEAD,fname=output_file, ns=ncols, nl=nrows, nb=1, interleave=0, data_type=4, offset=0, map_info=map_info, /write, /open
  ENDIF

  ;保存模型结果
  output_file = DIALOG_PICKFILE(TITLE = 'NDVI/npp与降水的偏相关系数T检验结果...')
  IF strlen(strtrim(output_file,2)) GT 0 THEN BEGIN
    ;二进制方式输出
    OPENW,lun, output_file,/get_lun
    ;写入npp结果数据到文件中
    WRITEU,lun,tchk_npppre
    FREE_LUN,lun
    ;写出头文件信息
    ENVI_SETUP_HEAD,fname=output_file, ns=ncols, nl=nrows, nb=1, interleave=0, data_type=4, offset=0, map_info=map_info, /write, /open
  ENDIF

  ;保存模型结果
  output_file = DIALOG_PICKFILE(TITLE = 'NDVI/npp与气温、降水的复相关系数F检验结果...')
  IF strlen(strtrim(output_file,2)) GT 0 THEN BEGIN
    ;ENVI_WRITE_ENVI_FILE, fchk_npptmppre, OUT_NAME = output_file,Map_Info = map_info
    ;二进制方式输出
    OPENW,lun, output_file,/get_lun
    ;写入npp结果数据到文件中
    WRITEU,lun,fchk_npptmppre
    FREE_LUN,lun
    ;写出头文件信息
    ENVI_SETUP_HEAD,fname=output_file, ns=ncols, nl=nrows, nb=1, interleave=0, data_type=4, offset=0, map_info=map_info, /write, /open
  ENDIF


  ;********************************结束****************************************
  ;print the end time
  endtime=systime(1)
  duration = double(endtime)- double(starttime)
  endstr = '-----持续时间:' + string(endtime)+ '-----'
  PRINT ,endstr

  Dlg_Result= Dialog_Message('相关分析执行完成(the analyst is finished)！',title='提示(tips)',/Information)
  ;销毁资源
  ;WIDGET_CONTROL, wgtlb, /destroy

  ;关闭envi二次开发模式
  ENVI_BATCH_EXIT
end







;***************************************************************
;discription： 复相关系数的F检验

;param1 = 列数
;param2 = 行数
;param3 = alpha,alpha=sqrt(n-m-1)
;param4 = 偏相关系数

;author: wangqiang cdut 2016-04-19
;email：  cqwensten@163.com  419436974@qq.com
;***************************************************************
FUNCTION func_check_f,ncols,nrows,alpha,frel_data
  ;f检验
  chk_result = make_array(ncols,nrows,type=4,value=0.0)

  print,'复相关系数-F检验计算中... ...',0
  steplen = float(nrows)/100

  FOR yindex=0,nrows-1 DO BEGIN
    FOR xindex=0,ncols-1 DO BEGIN
      frel_cellvalue = frel_data[xindex, yindex]
      IF(frel_cellvalue EQ -999.0) THEN BEGIN
        chk_result[xindex, yindex] = -999.0
      ENDIF ELSE BEGIN
        chk = alpha * ( frel_cellvalue^2.0 /sqrt(1.0 - frel_cellvalue^2.0 ))
        chk_result[xindex, yindex] = chk
      ENDELSE
    ENDFOR
    indexOfRunning = fix(yindex / steplen )
    print,'复相关系数-F检验计算中... ...',indexOfRunning
    WAIT,0.01
  ENDFOR

  print,'复相关系数-F检验计算中... ...',100
  return,chk_result
END


;***************************************************************
;discription： 偏相关系数的T检验

;param1=列数
;param2=行数
;param3=alpha,alpha=sqrt(n-m-1)
;param4=偏相关系数

;author: wangqiang cdut 2016-04-19
;email：  cqwensten@163.com  419436974@qq.com
;***************************************************************
FUNCTION func_check_t,ncols,nrows,alpha,prel_data
  ;t检验
  chk_result = make_array(ncols,nrows,type=4,value=0.0)

  print,'偏相关系数-T检验计算中... ...',0
  steplen = float(nrows)/100

  FOR yindex=0,nrows-1 DO BEGIN
    FOR xindex=0,ncols-1 DO BEGIN
      prel_cellvalue = prel_data[xindex, yindex]
      IF(prel_cellvalue EQ -999.0) THEN BEGIN
        chk_result[xindex, yindex] = -999.0
      ENDIF ELSE BEGIN
        chk = alpha*(prel_cellvalue /sqrt( 1.0 - prel_cellvalue^2.0 ) )
        chk_result[xindex, yindex] = chk
      ENDELSE
    ENDFOR
    indexOfRunning = fix(yindex / steplen )
    print,'偏相关系数-T检验计算中... ...',indexOfRunning
    WAIT,0.01
  ENDFOR
  print,'偏相关系数-T检验计算中... ...',100

  return,chk_result
END

;***************************************************************
;discription： npp与气温、降水的复相关系数

;author: wangqiang cdut 2016-04-19
;email：  cqwensten@163.com  419436974@qq.com
;***************************************************************
FUNCTION func_frel_b2b,ncols,nrows,npp_data,tmp_data,pre_data
  ;复相关系数
  frel_data = make_array(ncols,nrows,type=4,value=0.0)

  print,'复相关系数计算中... ...',0
  steplen = float(nrows)/100
  varible_arrary = make_array(ncols,nrows,2,type=4)
  FOR yindex=0,nrows-1 DO BEGIN
    FOR xindex=0,ncols-1 DO BEGIN
      npp_cellvalue = npp_data[xindex, yindex, 0]
      tmp_cellvalue = tmp_data[xindex, yindex, 0]
      pre_cellvalue = pre_data[xindex, yindex, 0]
      ;无数据区
      IF (npp_cellvalue EQ -0.9999)||(npp_cellvalue EQ -999.0)||(npp_cellvalue EQ 0.0)||(tmp_cellvalue EQ -999.0)||(pre_cellvalue EQ -999.0) THEN BEGIN
        frel_data[xindex, yindex] = -999.0
      ENDIF ELSE BEGIN
        npp_array = npp_data[xindex, yindex, *]
        tmp_array = tmp_data[xindex, yindex, *]
        pre_array = pre_data[xindex, yindex, *]

        npps = reform(npp_array)
        tmps = reform(tmp_array)
        pres = reform(pre_array)

        var_array = TRANSPOSE([[tmps],[pres]]) ;转置

        ;计算NPP与TMP的复相关性系数
        frel_npptmppre = M_CORRELATE(var_array,npps)
        frel_data[xindex, yindex] = frel_npptmppre
      ENDELSE
    ENDFOR

    indexOfRunning = fix(yindex / steplen )
    print,'复相关系数计算中... ...',indexOfRunning
    WAIT,0.01
  ENDFOR
  print,'复相关系数计算中... ...',100

  return, frel_data
END


;***************************************************************
;discription： npp与气温的偏相关系数（排除降水因子影响）

;author: wangqiang cdut 2016-04-19
;email：  cqwensten@163.com  419436974@qq.com
;***************************************************************
FUNCTION func_prel_b2b,ncols,nrows,npp_data,tmp_data,pre_data
  ;偏相关系数
  prel_data = make_array(ncols,nrows,type=4,value=0.0)

  print,'偏相关系数计算中... ...',0
  steplen = float(nrows)/100
  FOR yindex=0,nrows-1 DO BEGIN
    FOR xindex=0,ncols-1 DO BEGIN
      npp_cellvalue=npp_data[xindex, yindex, 0]
      tmp_cellvalue=tmp_data[xindex, yindex, 0]
      ;无数据区
      IF (npp_cellvalue EQ -0.9999)||(npp_cellvalue EQ -999.0) ||(npp_cellvalue EQ 0.0)|| (tmp_cellvalue EQ -999.0) THEN BEGIN
        prel_data[xindex, yindex] = -999.0
      ENDIF ELSE BEGIN
        npp_array = npp_data[xindex, yindex, *]
        tmp_array = tmp_data[xindex, yindex, *]
        pre_array = pre_data[xindex, yindex, *]

        ; Compute the partial correlation of X0 and X1 with the effects
        ; of X2 removed.
        ;result = P_CORRELATE(X0, X1, X2)
        ;计算NPP与TMP的偏相关性系数
        prel_npptmp = P_CORRELATE(reform( tmp_array),reform( npp_array),reform(pre_array))
        prel_data[xindex, yindex] = prel_npptmp
      ENDELSE
    ENDFOR
    indexOfRunning = fix(yindex / steplen )
    print,'偏相关系数计算中... ...',indexOfRunning
    WAIT,0.01
  ENDFOR
  print,'偏相关系数计算中... ...',100

  return,prel_data
END

;***************************************************************
;discription： npp与气温的Perison相关系数

;author: wangqiang cdut 2016-04-19
;email：  cqwensten@163.com  419436974@qq.com
;***************************************************************
FUNCTION func_relativity_b2b,ncols,nrows,npp_data,tmp_data
  ;Perison相关系数
  relativity = make_array(ncols,nrows,type=4,value=0.0)

  print,'线性相关系数计算中... ...',0
  steplen = float(nrows)/100
  FOR yindex=0,nrows-1 DO BEGIN
    FOR xindex=0,ncols-1 DO BEGIN
      npp_cellvalue=npp_data[xindex,yindex,0]
      tmp_cellvalue=tmp_data[xindex,yindex,0]
      ;无数据区
      IF (npp_cellvalue EQ -0.9999)||(npp_cellvalue EQ -999.0)||(npp_cellvalue EQ 0.0)||(tmp_cellvalue EQ -999.0) THEN BEGIN
        relativity[xindex,yindex] = -999.0
      ENDIF ELSE BEGIN
        npp_array =npp_data[xindex,yindex,*]
        tmp_array =tmp_data[xindex,yindex,*]
        ;计算NPP与TMP的线性相关性
        relativity[xindex,yindex] = CORRELATE(tmp_array,npp_array)
      ENDELSE
    ENDFOR

    indexOfRunning = fix(yindex / steplen )
    print,'线性相关系数计算中... ...',indexOfRunning
    WAIT,0.01
  ENDFOR
  print,'线性相关系数计算中... ...',100

  return,relativity
END
